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This paper extends our recent theoretical work concerning the feasibility of stable and accurate 
computation of turbulence using a large eddy simulation [Ida and Taniguchi, Phys. Rev. E 68, 
036705 (2003)]. In our previous paper, it was shown, based on a simple assumption regarding the 
instantaneous streamwise velocity, that the application of the Gaussian filter to the incompressible 
Navier-Stokes equations can result in the appearance of a numerically unstable term that can be 
decomposed into positive and negative viscosities. That result raises the question as to whether an 
accurate solution can be achieved by a numerically stable subgrid-scale model. In the present paper, 
based on assumptions regarding the statistically averaged velocity, we present similar theoretical 
investigations to show that in several situations, the shears appearing in the statistically averaged 
velocity field numerically destabilize the fluctuation components because of the derivation of a 
numerically unstable term that represents negative diffusion in a fixed direction. This finding can 
explain the problematic numerical instability that has been encountered in large eddy simulations 
of wall-bounded flows. The present result suggests that this numerical problem is universal in large 
eddy simulations, and that if there is no failure in modeling, the resulting subgrid-scale model can 
still have unstable characteristics; that is, the known instability problems of several existing subgrid- 
scale models are not something that one may remove simply by an artificial technique, but must be 
taken seriously so as to treat them accurately. 

PACS numbers: 47.11.-rj, 47.27.Eq, 47.10.+g, 83.85.Pt 
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I. INTRODUCTION 



Turbulence is one of the unsolved problems of physics 
0. Because a complete theoretical description has not 
yet been achieved even for a relatively simple flow con- 
figuration, numerical simulations are commonly used as 
a powerful tool for analyzing turbulent flows. Typi- 
cal numerical approaches include the direct numerical 
simulation (DNS), the Reynolds-averaged Navier-Stokes 
(RANS), and the large eddy simulation (LES). In DNS, 
all the scales of motion in the turbulent flows are re- 
solved by sufficiently fine computational grids, whereas 
in RANS, only the evolution of mean quantities is solved. 
LES is an intermediate technique between these ap- 
proaches, directly solving the large scales but modeling 
small-scale eddies by employing a subgrid-scale (SGS) 
model (or a subfilter-scale model) that approximately ac- 
counts for the effects of the small scales on the large scales 
0. Because LES enables us to solve time-dependent 
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large-scale turbulent flow problems with a relatively 
small computational time and storage compared to those 
required for DNS, this technique has recently been used 
not only for academic studies but also for practical in- 
dustrial flow computations that need time-dependent so- 
lutions. 

One of the major problems of LES is numerical insta- 
bility. As is already known, several existing SGS mod- 
els (e.g., the tensor-diffusivity 0, the scale-similarity 3, 
and the dynamic Smagorinsky p| models), which have 
been constructed based on a filtering procedure and the 
statistical properties of turbulence, have a numerically 
unstable property, and hence some artificial numerical 
treatments (e.g., smoothing or clipping of the SGS stress) 
have been incorporated so as to guarantee numerical sta- 
bility [E IE IE IE 111 • While the mechanisms of these mod- 
els' unstable behaviors have been described in the liter- 
ature (e.g., Refs. 0, IE IE El EI): the underlying reason 
as to why the SGS models have an unstable property 
has, to the authors' knowledge, not yet been fully clar- 
ified. To construct an excellent SGS model that is free 
from artificial, unphysical numerical treatments and has 
applicability to a wide range of flow configurations with 
high accuracy and robustness, it is necessary to elucidate 
the underlying mechanism of those unstable properties. 
Is this numerical problem caused by a failure in mod- 
eling or by other factors? To answer this question, it 
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should be meaningful to consider a similar but idealized 
question whether a completely accurate SGS model, if it 
exists, would be numerically stable. 

In Ref. ^} 1 Leonard has shown that the tensor- 
diffusivity model, which was derived by truncating an 
exact expansion series of the SGS stress terms and is 
thus exact under a certain condition, should behave un- 
stably along the stretching directions of fluid motion. 
This unstable behavior results from the negative diffusiv- 
ity of that model, which makes the governing equations 
ill-conditioned and leads to numerical instability when 
treated numerically by, e.g., finite difference methods. 
Winckelmans et al. @ have performed several numerical 
experiments using the pure (and mixed) tensor-diffusivity 
model and have pointed out that the model behaves un- 
stably under a wall-bounded flow condition, while it can 
provide a stable result for turbulent isotropic decay. In 
a recent paper 0] we have shown theoretically that the 
filtering procedure, which is the most fundamental com- 
ponent of LES, is a potential seed of the numerical in- 
stability in LES. Without having resorted to any SGS 
model, but based on a simple assumption regarding the 
streamwise component of flow velocity, we found that un- 
der a wall-bounded flow condition, the application of a 
Gaussian filter (one of the filters commonly used in LES 
to the Navier-Stokes equations results in the appear- 
ance of a numerically unstable term, a cross-derivative 
of the filtered velocity component, which can be decom- 
posed into diffusion and negative-diffusion terms. The 
above findings indicate that apparently the application 
of Gaussian filtering stabilizes the physical properties of 
the governing equations, but this is not always the case. 
That result implies that even a completely accurate SGS 
model that perfectly reproduces the physical properties 
of the true SGS components might not always be numer- 
ically stable. 

However, both Leonard's and Ida and Taniguchi's the- 
oretical works are insufficient to fully clarify the mech- 
anism of numerical instabilities in actual simulations. 
While those studies have only considered the instanta- 
neous nature of the SGS model or of the SGS forces, it 
is possible that their time-averaged nature is dissipative 
and that stable simulations can thus be achieved. In- 
deed, as mentioned above, the tensor-diffusivity model, 
which has an unstable property, canprovide stable solu- 
tions for an isotropic turbulent flow |8( ■ This observation 
suggests that further efforts must be made to improve 
our understanding of the numerical instabilities in LES. 

In the present paper, we have extended our previous 
work so as to obtain more acceptable and reliable results 
and to show that this numerical problem is universal in 
LES based on Gaussian filtering. We have previously 
assumed that the instantaneous value of the streamwise 
velocity component is linearly proportional to the dis- 
tance from a plane wall parallel to the bulk flow. In con- 
trast, in the present paper our discussion assumes that 
the statistically averaged streamwise velocity component 
is linearly proportional to the distance from the wall, an 



assumption that is more realistic because this is the case 
of a viscous sublayer forming near a plane wall. Fur- 
thermore, we discuss cases where an axisymmetric swirl 
exists in the statistically averaged velocity field. As will 
be shown in the following section, a numerically unstable 
term with a time-independent coefficient can appear in 
both situations. This term is always unstable in a fixed 
range of directions, while the previously discussed unsta- 
ble characteristics of the model and of the SGS terms 
depend on the directions of instantaneous stretching 
or of instantaneous shears [lj. Therefore, the present 
theoretical result can more accurately explain the nu- 
merical instabilities frequently encountered in inhomo- 
geneous flows that involve a strong steady shear in the 
mean velocity field. 

This paper is organized as follows: In Sec. Ill Bl we reex- 
amine the wall-bounded flow case, and in Sec. Ill Gl we ex- 
tend our theory to swirling flow cases. In the appendixes, 
we provide additional notes on further extensions of our 
theory. Section IlIII presents concluding remarks. 

II. THEORETICAL INVESTIGATIONS 

A. Filtering approach 

Incompressible, viscous fluid flows are described by the 
Navier-Stokes equations, which read: 

du t dUjUi dp d 2 u t . 

-ST + = + v for % = 1,2,3, 1 

Ot OXj OXi OXjOXj 



where Einstein's summation convention is assumed, and 
Ui = Ui(x\,X2,xz,t) are the velocity components, p is 
the pressure divided by the constant fluid density, and v 
is the kinematic viscosity. In LES, a filter is applied to 
this system of equations to separate the large and small 
scales. This filtering procedure is in general achieved by 
the following convolution: 

X=oo 

F{x,...,t)= J L(x-X)F(X, ...,t)dX, (3) 

where (•) denotes the filtered quantity, x is an indepen- 
dent variable of an arbitrary function F, and L{X) is the 
filter function. In the present study, we assume L(X) to 
be the Gaussian function 

which satisfies fx=™ L(X)dX = 1, where A is the filter 
width (assumed to be constant), and 7 is a real, positive 
constant. In the previous paper we set to 7 = 1/2 as 
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done by Klimas for the Vlasov equation, while 7 = 6 
is generally used in LES (l3|: we adopt the latter value 
in the present study. Appling the filtering operation to 
Eqs. (HJ and © yields 



du t 
dt 



dujUi 
dxi 



dp 



d 2 u t 



dxi dxjdxj 



(5) 



duj 
dt 



dujUi 
dxj 



dp 

dxi 



d 2 Ui 



dxj dxj dxj 



^ (6) 



with 



where we used 



duj 
dxi 



= 0, 



(7) 



df\_df fdf\ df 



dt 



Of 



dx„ 



dxj 



(/ being a dependent variable), and r»j is the so-called 
SGS stress tensor that generally needs to be modeled. 

Before getting into the main subject, we present here 
some mathematical tools useful for the present investiga- 
tion. Supposing that the Gaussian filter is applied in the 
Xi direction, we have 



(Xif) =Xif + 



a 2 af_ 

27 dxi 



and 



(xjf)=Xjf forj ¥= h 



(8) 



(9) 



which can be g iven in the same manner as shown in 
Refs. [l2lll^ . lT5| . and have been utilized in, e.g., Ref. [It| . 
Using Eq. (JSJ, we have 



Xi Xi 



(*?) 



A 2 

27 • 



(10) 



(11) 



Based on these results, in the following subsections we 
perform theoretical investigations concerning the Gaus- 
sian filtered Navier-Stokes equations under a plane chan- 
nel and swirling flow conditions. 



B. Plane channel flow 

Suppose that x\, X2, and 23, respectively, are the 
streamwise, the wall-normal, and the spanwise directions, 
and a plane solid wall is set for X2 < 0. Decomposing the 



velocity components into the statistically averaged value 
Ui and fluctuation part u^, we have 



Ui{x,t) = Ui(x) +w-(x,i), 



(12) 



where Ui is assumed to not depend on time. Further- 
more, we assume that 



U(x) - (C/ 1 (x),C/ 2 (x),[/ 3 (x)) 
= {(3x 2 ,0,0) 



(13) 



with /3 being a real constant, i.e., the mean streamwise 
velocity is linearly proportional to the distance from the 
wall. (The extension to more general cases is briefly dis- 
cussed in Appendix 1X1) Based on these assumptions, one 
knows 



dUj du'j 
dxi dx, 



= 0. 



(14) 



Substituting Eqs. (|12fl and (|13fl into the convection 
terms in the equation for u\ and using Eq. I|14|l , we have 



dujUi _ du{ dUi , du'ju'-t 

U\— 1 Un ' 



dxi 



dxi dx2 



dxj 



= (3x2^ + 13^2 + ^—. 
dxi dxj 



(15) 



Appling Gaussian filtering in the wall-normal direction 
to this equation yields 



du]W - du[ A 2 d 2 u' dUi , du'M 



dxi 



dxi 27 dx\dx2 dx2 



dxi 



du[ A 2 d 2 u[ „ , du' jU [ 
(3x2^ + f3— — ^ + I3u' 2 1 3 



' dx\ 27 dx\dx2 



dxj 



(16) 



where we use Eq. and U\ = U\ given easily by 
Eq. (|10|l . Even if a three-dimensional Gaussian filter is 
adopted, almost the same formula is derived because U\ 
only depends on X2, and hence no additional term is de- 
rived by filtering in the other directions. Meanwhile, the 
convection terms in terms of the filtered velocity compo- 
nents, d{ujU\)/dxj, can be written as 



dujUx - du[ dUx _, du'^ 
— Ui— — 1 ' 



dxj 



dx\ dx2 



dxj 



(17) 



From Eqs. (|16fl and (|17|) . we have 



d nj = p A 2 d 2 u[ + d^v{-v^) (ig) 



dxj 27 dx\dx2 



dxj 



Substituting this and d 2 U\/ dxjdxj = into Eq. © 
yields 



dui dujUi dp d 2 u[ A 2 d 2 u[ 



dt dxj 



0- 



dx\ dxjdxj 27 dx\dx2 



dtyui - UjM'i) 
dx i 



(19) 
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which can be considered an equation for both Ui and u[ 
because 

dui du' x 
~dT = ~W' 

The second-to-the-last term of Eq. (I19fl is a cross- 
derivative of the dependent variable u'-, . As has been 
proved in several studies (e.g., Refs. [H,Il5j]), the deriva- 
tives of this type should be numerically unstable when, 
for instance, a finite difference technique is adopted for 
solving the equation. In what follows, we consider the 
stability condition for Eq. (|19|l . 

The last three terms of Eq. 1|19|1 may make a dominant 
contribution to the numerical stability of LES. We will 
first focus our attention on the second- and third-to-the- 
last terms: 



d 2 u[ 

dxjdxj 



A 



d 2 u\ 



27 dx\dx2 

Appling the coordinate transformation of 



%\ + X2 —X\ + Xi 



V2 



V2 



(20) 



(21) 



which represents rotation by 45° around the X3 axis, 
Eq. H2U|) is rewritten as 



<9£ 2 drf 



Oxl 



-13 



A 2 



2 7 V 2 dSf 2 drf 



,(22) 

From this, it can be seen that the following restriction is 
required to ensure that Eq. ll'L'l) is a positive viscosity: 



A 2 



(23) 



The shear gradient (3 may have different values de- 
pending on the problems considered. We consider here a 
simple example to show concretely how restriction <|23[) 
works in a realistic situation. As is well known, in the vis- 
cous sublayer of turbulent channel flo ws, the statistically 
averaged streamwise velocity follows 17j 



/ \ / 14-7- 

Ui{x 2 ) = u T x 2 

v 



(24) 



with u T being the wall-friction velocity. Using Eq. 124|) 
and rewriting reduce Eq. ill' Ml) to 



—A = A + < 2^, 



(25) 



where the superscript + denotes a quantity in the wall 
units. For 7 = 6, this result essentially corresponds to the 
stability condition derived by Kobayashi and Shimomura 
from the tensor-diffusivity term in a dynamic SGS model 

We should note here that Eq. I|25|) is a minimum condi- 
tion for stabilizing the numerical solution; the last term 



of Eq. Q19[l . which we have neglected, can also be a seed 
of numerical instability because if, for example, oc x 2 
is true at a certain instant and location, then 



1 dx 1 



_, du[ d 2 u[ 
1 dxi dx\dx 2 



(See Ref. Nevertheless, this incomplete condition 

creates a strong restriction on the grid width (h). The 
grid width should satisfy h+ < ~ 2.4 for A+ = 2h + , 
a setting that has frequently been used, or h + < V6/2 ~ 
1.2 for A + = Ah + , which is necessary for the contribution 
of the SGS force to be significant compared to the trunca- 
tion and aliasing errors of a second-order finite difference 
scheme [I3 ■ (We should note here that this restriction is 
imposed only on the grid width in the wall-normal direc- 
tion because only the filtering in this direction yields the 
cross-derivative term, as mentioned already.) This result 
implies that at least several grid points are needed in the 
viscous sublayer to guarantee numerical stability, which 
is a formidable restriction in practical applications of a 
large-scale high-Reynolds-number turbulent flow. It is 
known that when the inner layer of turbulent boundary- 
layer flows is resolved, the number of grid points required 
for LES exceeds current computational capacities even at 
moderate Reynolds numbers, and hence developments of 
SGS models that are applicable to a coarse wall-normal 
grid and of wall models that do not need a strong refine- 
ment of the near-wall grids are the most pressing issues 
of LES; see a recent review by Piomelli and Balaras [l9j . 
The numerical problem that we have suggested above 
provides a taste of the difficulties of constructing an ex- 
cellent model. 

The stability condition presented above would be 
partly weakened if the last term of Eq. I|19|) is mod- 
eled using an eddy-viscosity model. However, this treat- 
ment does not always achieve a stable solution. Actually, 
Winckclmans ct al. [8| have observed by numerical exper- 
iments that for a plane channel flow, a dynamic mixed 
model based on the tensor-diffusivity model (which can 
accurately represent the characteristics of second-order 
cross derivatives [TT||') generates numerical instabilities 
that eventually make the numerical simulation blow up. 
To remove those instabilities, they were forced to add 
an artificial damping coefficient to the tensor-diffusivity 
portion. In the mixed model that they used, the coeffi- 
cient of eddy-viscosity part was well tuned by a dynamic 
procedure. Essentially the same observations of the nu- 
merical instability of the dynamic mixed model can be 
found in Ref. |ll| by Kobayashi and Shimomura. In those 
studies, the wall-normal grid width was set to be sev- 
eral times larger than that used in DNS. That is, those 
studies suggest that, for a large grid (and filter) width, 
which is much-needed in practical applications, the ad- 
dition of an eddy-viscosity term does not always guar- 
antee numerical stability. If, on the other hand, all the 
SGS terms are modeled using an eddy-viscosity model 
whose coefficient is guaranteed to be positive, the nu- 
merical stability of the filtered system would be guaran- 
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teed in many cases. However, it has been suggested by 
many researchers that eddy-viscosity models have several 
deficiencies that make the LES results suspicious (e.g., 

Refs. dHHlillHlililliil). For instance, as 
has frequently been demonstrated by a priori tests us- 
ing DNS techniques, the SGS quantities modeled using 
eddy-viscosity models do not correlate well with the ac- 
tual SGS quantities |2j, 111 |2Jj, |23|, |26( , implying that these 
models have low accuracy and/or narrow applicability. 
Indeed, in the present case the eddy-viscosity models 
cannot describe the cross-derivative term in Eq. (|19fl . 
which is a hybrid of positive and negative diffusions. 
Moreover, eddy- viscosity models were constructed under 
the assumptions that the subgrid turbulence is isotropic 
and the filter scale is in an inertial range, assump- 
tions that are in general both violated in inhomoge- 
neous turbulent flows, especially near walls. Further- 
more, the dynamic types of eddy-viscosity models need 
artificial techniques for smoothing the model coefficient 
and for clipping its negative values 0, IE IE E3- 
Based on the above, we consider the eddy-viscosity mod- 
els to be insufficient. Most of the recent efforts at 
SGS modeling have been devoted to developing a model 
that does not stron gly rely on eddy- viscosity models 
0, IE HE ISE EE 011 EL Models based on a kind 
of "defiltering" procedure |25l I29I l30| are one such ap- 
proach, determining significant portions of the SGS terms 
by an analytical, not empirical, procedure. The mixed 
models based on the tensor-diffusivity term 0, IE Hill can 
also be categorized into this type |2g. We should note 
here that attempts to improve eddy-viscosity models are 
also being continued; see, e.g., Refs. 0i EE HE H3 ■ 

It is worth noting that the numerical stability of the 
cross-derivative term in Eq. I|19[l is time-independent, as 
this term always reveals a negative diffusivity in a fixed 
range of directions, whereas the numerical stability of the 
last term of the same equation may be time-dependent. 
As mentioned in Sec. [I] Winckelmans et al. have pointed 
out that the numerical stability of the tensor-diffusivity 
model depends on the problems that they have consid- 
ered; obtaining stable results is difficult to achieve for a 
strongly inhomogeneous (wall-bounded) flow, even when 
an eddy-viscosity term is added, but is possible for a ho- 
mogeneous flow even by the pure tensor-diffusivity model 
Q. They have explained these results as follows: For 
the homogeneous flow, the direction in which the tensor- 
diffusivity term reveals the negative diffusivity evolves 
continuously in time and space, and hence the model 
term may be dissipative on time average. However, in 
the inhomogeneous flow case, long-lived negative diffu- 
sion events can take place near the wall, leading to nu- 
merical instability and a resulting divergence of the solu- 
tion. As mentioned above, those numerical findings indi- 
cate an inadequacy of the previous stability analyses by 
Leonard ^(J an d by Ida and Taniguchi which have 
only considered the instantaneous behavior of the SGS 
model or the SGS term. The present theoretical find- 
ings, showing the appearance of a term that is always 



numerically unstable in fixed directions and hence is un- 
stable on time-average as well, can more correctly explain 
the problematic numerical instability encountered in the 
simulations of plane channel flows [E El ■ 



C. Swirling flow 

Turbulent flows involving a large-scale swirl are of 
practical importance in connection with combustion en- 
gineering, aeroacoustics, meteorology, and so on. In this 
subsection, we discuss the LES of turbulent flows having 
a rotating mean velocity about the x% (or z) axis, and 
show that the existence of the swirl can cause numerical 
instability through the Gaussian filtering operation. 

First, we consider a very simple case in which the swirl 
has a constant angular velocity; that is, the mean velocity 
takes the form of 

U(x) = (oara.-OKBi.O), (26) 

where a is a real constant. This mean velocity satisfies 

U • x* = 0, 



|U|/|x*| = |a|, 



and 



dUj 
dx t 



0, 



where x* = (x\, X2, 0). As a result, we know 
du' 



dx, 



= 0. 



(27) 



Using the above, the convection terms for u% are rewrit- 
ten as 



dxj 



= U 1 — — - + [T 2 — -i + —tvL + U 2 ^ r - L + 



dxi 8x2 dx 



8x2 dxj 



(28) 

where the first four terms on the right-hand side (RHS) 
result from the existence of the swirl, among which the 
first and second terms, 

du\ du\ du\ du\ ,. 

Ui-± + U 2l r±=ax2 1 r±-ax l --±, (29) 
oxi ox 2 0x1 0x2 

may have a significant influence on the numerical stability 
of the equation for Uj . 

We show here that the filtered formula and its stabil- 
ity change depending on how the Gaussian filtering is 
applied. Appling the Gaussian filter in the X2 direction 
to Eq. results in 



0x1 J \ 0x2 J 0x1 2-f 0x10x2 0x2 
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which can be rewritten as 



Ui 



du[ 



du[ _ ^A 2 d 2 u[ 

3X2 



27 dx\dx2 
(30) 

The cross-derivative term appearing on the RHS may 
cause numerical instability when the absolute value of 
its coefficient is sufficiently large. In contrast, when the 
Gaussian filter is applied in both the x\ and X2 directions, 
we obtain a stable formula, 



du[ 
dx\ 



U 



dx 2 



ax 2 



-ax\ 



ax2 



du[ 
dx\ 
du[ 
8x2 
du[ _ 
dx\ 



a 



d 2 u[ 



27 dx\dx2 
A 2 d 2 u[ 



a 



ax 1 



27 dx\dx2 
du[ 



8x2 



and thus, 



Ui 



dxi 



OX\ 



u 



du[ 

3X2 



u 



dx 2 



0. (31) 



Equation l|31|) has no cross derivative. 

Next, we consider the same example but in the cylin- 
drical coordinate (r, 0, z). In this case, the mean velocity 
(|26[1 and the instantaneous velocity, respectively, are rep- 
resented by 



and 



U(r, 0, z) = (U r , U B , U z ) = (0, -ar, 0), (32) 



(u r ,u g ,u z ) = (u' r ,Ug + u e ,u' z ), (33) 



and the nonlinear terms in the Navier-Stokes equation 
for ug take the form of 



u r - 



dug 



Ug dug 

7¥ 



du 



e , u r ug 



dr r dO dz r 

where the velocity components should satisfy 

u r 1 dug 
~ + r~dT 



(34) 







8u r 

dr 
du' r 
dr 



r 



du z 
dz 
du' 



I du' g 

r d6 dz 



Substituting Eq. lO into Eq. ((331 yields 
Ug du'g 



du 'g 
dr 



[dUg 


, + Ue 


\ dr 


r 


+ ^± 


dUg 


r 


~de + 



./ du'g 

z dz 



u„u e 



(36) 



The terms in the last parentheses represent the nonlinear 
interaction between the fluctuation components, whereas 
the remaining terms describe the interaction between the 
mean and fluctuation portions. Appling the Gaussian 



filter in the r direction with a filter width of A r 
first term results in 



Ug du'g 
r dO 



du'g 



Ug du'g 
r dO 



to the 



(37) 



where we used Eq. I|10|l and assumed r ^> A r so that the 
value of the Gaussian filter is sufficiently close to zero 
at r = 0. No unstable term appears in Eq. (|37|l . Even 
when the Gaussian filter is applied in both the r and 
directions, almost the same formula is obtained because 
Ug depends only on r. Specifically, in the present case 
the resulting formula is not, unlike the previous case, 
dependent on how the Gaussian filter is applied. 

The above results, which reveal that the numerical 
stability of LES depends not only on the flow configu- 
ration but also on the filtering strategy and coordinate 
system adopted, can be interpreted as follows: As has 
been shown in Ref. 01 an d m the preceding subsection, 
the existence of a shear leads to the appearance of a nu- 
merically unstable term. Although no shear exists in the 
mean flow profile described by Eq. I|2t)|) , a virtual shear is 
observed when the Gaussian filter is applied only in one 
direction in the Cartesian coordinate (x\ or X2), resulting 
in the derivation of the cross-derivative term. This prob- 
lem does not arise for the cylindrical coordinate, since 
the mean convection velocity in the (r, 0) space, 



dr dO 
~dt' ~dt 



U r ^) =((>.-„) 
r 1 



is uniform. 

Lastly, we consider a case where the mean velocity has 
a second-order term, i.e., assuming 



U(r,6»,z) = (0,air + a 2 r 2 ,0), 



(38) 



where ot\ and 0,2 are real constants. Substituting this 
into the first term on the RHS of Eq. (|36|l and applying 
the Gaussian filter in the r direction, we have 



Ug du'\ , , du'g A 2 d 2 u' , . 



which can be rewritten as 



(35) (Ue&u'e} Ue du' g _ ^ A£c% 
d0 



Ug dvJg 

d0 



_ A 2 du' g a; o"u' e 
a2 ¥^^ +a2 2^dT Z dO 



(40) 



because, based on Eqs. (|rU|) and l(TT|) . 



Ug = air + a 2 [ r 2 + — - 



A 2 

27 

The last term of Eq. H4U|) is a cross derivative of the 
dependent variable and should destabilize the numerical 
solution, whereas the second-to-the-last term represents 
convection in the direction that can be solved stably. 
For a2 = 0, the RHS terms of Eq. (|4*U)) vanish, and this 
equation corresponds to Eq. I|37|l . More specifically, the 
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instability in this case is due to the second-order compo- 
nent in the mean velocity, which violates the uniformity 
of the angular velocity and thus represents a shear. 

As in the plane-channel case, the cross-derivative terms 
derived above have a timc-indcpcndcnt coefficient and 
thus should lead to numerical instability if an unsuitably 
large filter is adopted. The present results imply that 
such an unstable term can appear in many situations. In 
a future paper, we will consider more general and compli- 
cated flow configurations to create a generalized theory 
of numerical instability due to filtering operations [23] ■ 



III. SUMMARY AND CONCLUSION 

In summary, we have theoretically investigated the nu- 
merical instability of LES caused by the filtering proce- 
dure based on simple assumptions regarding the statis- 
tically averaged velocity, and have shown that Gaussian 
filtering yields a numerically unstable term, which always 
reveals a negative diffusivity in a fixed range of directions 
in cases of both plane channel flow and swirling flow. 
This conclusion confirms and extends that of our previ- 
ous work ^3] • We anticipate that it would be interesting 
to examine the relationship between this numerical in- 
stability and the shear-induced (physical) instabilities in 
turbulence [Sfillsfil]. 

Furthermore, we have presented several additional 
findings. In Sec. Ill Bl we showed that in channel- 
flow case, this numerical problem strongly restricts the 
wall-normal grid width necessary for achieving stable 
and accurate solutions. In the academic computa- 
tions of a plane channel flow using an LES technique, 
the grid width in the wall-normal direction has some- 
times been set to almost the same as that required in 
DNS, and no filter is adopted in this direction (e.g., 
Refs. H HJ El HI the present result provides a 

grounding for this custom. Furthermore, this result in- 
dicates a significant difficulty in developing an excellent 
SGS model that can provide satisfactory (i.e., not only 
stable but also accurate) solutions, even with large wall- 
normal filter widths, which is seriously desired in practi- 
cal applications of high- (and even moderate-) Reynolds- 
number turbulence. In Sec. Ill CI based on investigations 
of swirling flow cases, we pointed out that the numeri- 
cal stability of the filtered equations depends not only on 
the flow configuration, but also on the dimension of the 
Gaussian filtering and the adopted coordinate system. 
This finding would be of particular importance in practi- 
cal applications that use a general curvilinear coordinate 
system or unstructured grids. 

In the present study (and also in our previous paper), 
it was assumed implicitly that all the spectral modes con- 
tained in the velocity field are fully resolved because the 
Gaussian filter does not cut off any modes, but only 
damps the high-frequency modes. (See Appendix [5] 
which presents a related remark concerning cases where 
the spectral cutoff or the top-hat filter is employed, both 



of which eliminate certain modes.) The present results 
therefore indicate that even in such an ideal case, the 
complete SGS model (or, more precisely, the complete 
subfilter-scale model whose characteristic length is de- 
termined independently from the grid width), which can 
completely reproduce the properties of real SGS stresses, 
can be numerically unstable when treated numerically, 
since the complete model perfectly describes the char- 
acteristics of the unstable cross-derivative term. This 
problem poses a dilemma for practitioners of LES who 
are looking for an accurate and stable SGS model. To 
improve model accuracy, this numerical instability prob- 
lem should be confronted. However, if a rough, artificial 
treatment is adopted for stability, the accuracy of the 
LES results will not be guaranteed. A potential approach 
to overcoming this difficulty is, as suggested in our pre- 
vious paper |12j . to construct a stable and accurate nu- 
merical solver for the numerically unstable term, though 
this would admittedly be quite a difficult task. (Leonard 
[To| and Moeleker and Leonard [37J have proposed La- 
grangian methods based on the tensor-diffusivity model 
and an anisotropic particle scheme to solve the negative- 
diffusion problems, and have achieved stable solutions for 
a two-dimensional scalar transport equation with known 
velocity fields. However, several issues remain to be over- 
come to extend those methods to the case of Navier- 
Stokes turbulence [37J ■) Further careful and vigorous 
considerations are necessary for this instability problem 
to be solved. 
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APPENDIX A: CLOSURE FOR HIGH-ORDER 
MEAN VELOCITY 

As a first step toward constructing a general theory for 
the filtering instability, let us consider the closure of 



dxi 



in the case where U\ is a general function of X2, under 
the assumption that the Gaussian filter is applied in the 
x-i direction. Using the Taylor expansion, the arbitrary 
function U\{x-i) is represented by a polynomial of infinite 
order: 



n=0 
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where a n are real constants. Therefore, what we have to 
do is to consider the closure of 



(Al) 



with n being an arbitrary positive integer. This aim is 
achieved as follows: Using Eq. JSJ successively, Eq. (|A1|) 



is rewritten into a closed form in terms of 



u\ , as 



X 2 x. 



1 ® u 'i 

1 di^ 



X\ x 



i-X du 'i 
' 2 d Xl 



where 



X 



xx yx 2 

X n : 

ax\ 



A 2 d 
27 dx 2 ' 



-2 9u{ 
dx\ 



(A2) 



The resulting formula has derivatives of high orders, and 
hence stability analysis would be somewhat intricate. 



APPENDIX B: SPECTRAL CUTOFF AND 
TOP-HAT FILTERS 

In Ref . [l4( , Klimas showed that the Gaussian filtered 
Vlasov equation can be rewritten into a closed form in 
terms of the filtered distribution function without any 
approximation. Also, in Ref. [3j| ( see also Refs. [ToL 12^1 
Yeo showed that the Gaussian filtered Navier-Stokes 
equations are written in a closed form having an infi- 
nite series. Their results have been utilized in our study, 
as described in the present and previous papers. This 
Appendix is devoted to showing that when the spectral 
cutoff or the top-hat filter is adopted, it is in general im- 
possible to analytically derive an exact, closed formula. 

Let us consider two arbitrary velocity fields 



(m) 



(4 m) (x) 



,( m ) 



(x),4 m) (x)) 



for m = 1, 2. 



whose spectra are 

u(" l '(k) = (4 m) (k),4 m) (k),4 m) ( k )) for m = X > 2 ' 

where x = {xi,x%,xa) and k is the wavenumber vector. 
If these velocities, at a certain instant t\, satisfy 



fiW(k) ^ u< 2 >(k) for |k| > k c 



but 



and 



•,(D 



(k) 



a (2) 



(k) for |k| < k c 



U W(x) =u( 2 )(x) for all x, 



where (•) denotes the spectral cutoff with an identical 
cutoff wavenumber k c , then the closed equation will pro- 
vide the same result after this instant [i.e., uW(x, t) — 

u( 2 )(x, t) for t > tt], although the unfiltered velocities 
have different high-wavenumber modes. This unaccept- 
able conclusion has resulted from the assumption that a 
closed formula exists. 
If, on the other hand, 



«(x)=i.( 2 > 



ir 



'(x) + cos(fcx 2 )i 



is true at an instant (where A; is a constant wavenumber 
and i denotes the unit vector in the x\ direction) and 
the wavelength of the cosine function in this equation is, 
for example, equal to the filter width of the top-hat filter 
applied (at least) in the x 2 direction, then 

uW(x) = ^(x) for all x. 

(Here, (•) denotes the top-hat filtering.) As in the pre- 
vious case, this result denies the existence of a closed 
formula. 

The above results imply that in order for an exact 
closed formula to exist, the spectrum of the applied filter 
needs to have a nonzero value for |k| 7^ 00. The convolu- 
tion product in the physical space using a filter function 
L, 



u 



M(x) 



L(x-Ou< m >(0de, 



is represented in the Fourier space by the simple multi- 
plication of the spectra of L and u' m ^: 



uW(k) = L(k)u( m '(k). 

Based on the above, we know that if uW(k) =^ u( 2 '(k) 
for certain k and L(k) has a nonzero value for all k, then 
fiW(k) 7^ u( 2 )(k) at least for the certain k and conse- 
quently uW(x) 7^ u( 2 '(x) at least in some regions; more 
specifically, if the adopted filter function is not equal to 
zero for all k, filtering different velocity fields results in 
different filtered velocity fields. 

The present conclusion is in opposition to the theo- 
retical result of Carati et al. |2^] which proved mathe- 
matically the existence of a closed formula in the top-hat 
case. Let us consider the cause of this discrepancy. In 
the study of Carati et al., it was first assumed that the 
generalized expansion series for one-dimensional filtering 
in the x direction takes the form of 



(ab) 



Crsdl.adib, 



(Bl) 



where a and b are arbitrary, continuous, and differen- 
tiable functions of x and c rs are real constants. We show 
here that this assumption is not always valid for top-hat 
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filters. If a is a sinusoidal wave whose wavelength is equal 
to the characteristic width of the applied top-hat filter, 
then a(x) = for any value of x and consequently the 
right-hand side of Eq. (|B1|1 is also equal to zero for any 
x if finite values. However, the left-hand side of 

this equation is not necessarily zero under the present 
condition; i f, for example, b = a, then ab > for any 
x and thus (ab) > 0. This contradiction is caused by 
the fact that the value of the top-hat filters can be zero 
in the Fourier space. The generating function, Eq. (3.7) 



of Ref. [26J used to derive the generalized expansion se- 
ries, is definable only for the filters that always have a 
nonzero value in the Fourier space as Gaussian filters, 
because that function diverges at the wavenumbers for 
which the filter value is zero. Therefore, the expansion 
series for top-hat filters is definable only if the filter width 
is smaller than the smallest resolved scale. When this re- 
striction is not fulfilled, the expansion series should not 
converge. 
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